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ABSTRACT 


A  steady,  axially  symmetric,  incompressible,  viscous  flow  past  a  rigid  sphere  is 
numerically  simulated  by  using  a  numerical  scheme,  based  on  spectral  methods.  The 
equations  have  been  reduced  to  two  sets  of  nonlinear  second  order  partial  differential 
equations  in  terms  of  vorticity  and  stream  function.  The  calculations  have  been  carried  out 
for  Reynolds  numbers,  based  on  the  sphere  diameter,  in  the  range  0.1  to  104. 

The  numerical  results  have  verified  that  there  is  excellent  agreement  with  Stokes 
theory  at  very  low  Reynolds  numbers.  At  moderate  to  intermediate  Reynolds  numbers 
there  is  good  general  agreement  with  available  experimental  data  and  flow  visualization 
pictures.  The  Reynolds  number  at  which  separation  occurs  is  estimated  as  20.  The 
approach  to  boundary-layer  behavior  with  increasing  Reynolds  numbers  is  also  verified  by 
comparison  with  potential  flow  theory  and  analytical  boundary-layer  solution. 
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I.  INTRODUCTION 


Numerical  simulations  in  fluid  dynamics  have  historically  attracted  considerable 
attention,  through  the  use  of  various  numerical  schemes  such  as  finite-difference  and 
finite-element  methods.  In  this  study,  the  direct  numerical  simulation  of  a  time-dependent, 
axially  symmetric,  incompressible  viscous  flow  past  a  rigid  sphere  was  studied  by  using  a 
different  numerical  scheme  based  on  spectral  methods.  Spectral  methods  have  become 
increasingly  popular  in  recent  years,  especially  with  the  development  of  fast  transform 
methods.  It  has  been  successfully  applied  to  numerical  weather  predictions,  numerical 
simulation  of  laminar/turbulent  flows  and  other  fundamental  problems  especially  where 
high  accuracy  is  desired  for  relatively  lower  memory  requirements. 

The  drag  dependence  on  Reynolds  number  was  studied  and  the  results  were 
compared  with  the  previous  results.  The  stream  function  \j/  and  the  vorticity  co  were 
visualized  through  the  use  of  the  numerical  results  and  compared  with  the  flow 
visualization  pictures  for  various  Reynolds  numbers.  The  boundary-layer  thickness  on  the 
surface  of  the  sphere  was  also  investigated  through  the  use  of  the  results  of  the 
simulation.  The  boundary-layer  thickness  on  the  surface  of  the  sphere  is  calculated  by 
introducing  the  power  series  expansions  of  the  stream  function  into  the  boundary-layer 
equations.  The  results  are  valid  for  Re  — >  <*>  analytically.  The  minimum  Reynolds  number 
for  which  the  boundary-layer  assumption  is  valid,  is  investigated  by  comparing  the 
boundary-layer  thickness  which  is  obtained  from  the  numerical  solution,  with  analytical 
results  in  the  boundary-layer  limit 

It  is  clear  from  the  foregoing  that  the  objective  of  this  study  is  to  carry  out  the 
numerical  scheme  based  on  the  spectral  methods,  through  the  use  of  the  vorticity-stream 
function  form  of  the  Navier-Stokes  equations,  on  a  time-dependent,  axisymmetric,  viscous 
flow  about  a  rigid  sphere.  A  Matlab  program  was  developed  to  implement  the  foregoing 
numerical  scheme.  The  results  were  compared  to  the  previous  studies  which  were  based 
on  different  methods.  However  the  details  of  the  method  itself  were  not  compared  to  the 


other  methods  which  were  used  in  previous  studies.  The  expectations  are  that  the  results 
will  match  with  the  previous  results,  improve  our  current  ability  to  compute/predict  the 
evolution  of  flow  for  a  particular  geometry  and  conditions,  and,  hopefully  shed  some  light 
on  the  physics  of  flows  which  have  been  heretofore  uncalculated. 


II.  BACKGROUND  STUDIES 


The  collocation  method  appears  to  have  been  first  used  by  Slater  (1934)  and 
Kantorovic  (1934)  in  specific  applications.  It  was  developed  as  a  general  method  for 
solving  ordinary  differential  equations.  In  1938  Lanczos'  work  showed  for  the  first  time 
that  a  proper  choice  of  trial  functions  and  distribution  of  collocation  points  is  crucial  to  the 
accuracy  of  solution.  This  method  was  revived  by  Clenshaw  (1957),  Norton  (1963)  and 
Wright  (1964).  The  application  of  Chebyshev  polynomial  expansions  to  the  initial  value 
problems  are  involved  in  these  studies. 

The  spectral  collocation  method  was  applied  to  partial  differential  equations  for 
spatially  periodic  problems  by  Kreiss  and  Oliger  ( 1 972)  (who  called  it  the  Fourier  method) 
and  Orszag  (1972)  (who  termed  it  pseudospectral).  These  were  the  earliest  applications  of 
spectral  collocation  or  pseudospectral  method.  This  approach  was  very  attractive  because 
of  its  application  to  variable-coefficient  and  even  non-linear  problems. 

The  Galerkin  approach  which  depends  on  the  same  trial  and  test  functions,  was 
applied  to  a  meteorological  model  by  Silberman  (1954).  This  was  the  first  serious 
application  of  spectral  methods  to  PDE's.  After  Orszag  (1969,1970)  and  Eliasen, 
Machenhauer  and  Rasmussen  (1970)  developed  transform  methods  for  evaluating 
convolution  sums  arising  from  quadratic  non-linearity,  spectral  Galerkin  methods  only 
became  practical  for  high  resolution  calculations  of  such  non-linear  problems. 

The  first  unifying  mathematical  assessment  of  the  theory  of  spectral  methods  was 
covered  in  the  monograph  by  Gottlieb  and  Orszag  (1977).  Since  then,  the  theory  has  been 
extended  to  cover  many  different  problems,  such  as  variable-coefficient  and  non-linear 
problems.  Applications  in  fluid  dynamics  were  reviewed  in  the  symposium  proceedings 
edited  by  Voigt,  Gottlieb  and  Hussaini  (1984).  Other  introductory  or  review  articles  have 
appeared  recently  such  as  Mercier  (1981),  Hussaini,  Salas  and  Zang  (1985).  Zang  and 
Hussaini  (1985),  Deville  (1984),  Gottlieb  (1985)  and  Hussaini  and  Zang  (1987). 


It  is  well  known  that  the  background  of  this  particular  problem  of  flow  over  a 
sphere  has  a  long  history  due  to  the  fundamental  interest  in  external  flow  over  bluff 
bodies.  The  first  theoretical  results  for  steady,  uniform  flow  past  a  rigid  sphere  were  given 
by  Proudman  and  Pearson  (1957)  and  are  valid  for  Re  <  1.  For  higher  Reynolds  numbers 
computations  have  been  conducted  by  Rimon  and  Cheng  (1969),  Le  Clair,  Hamielec  and 
Pruppacher  (1970)  and  more  recently  by  Fornberg  (1988).  All  these  computations  are 
based  on  finite-difference  schemes.  Dennis  and  Walker  (1971),  Oliver  and  Chung  (1985) 
and  Brabston  and  Keller  (1975)  used  a  different  computational  approach  based  on  an 
expansion  of  Legendre  functions  and  spherical  harmonics.  The  results  of  these  simulations 
are  consistent  with  the  experimental  data  for  the  drag  forces,  and  have  been  used  to 
develop  improved  formulae.  Marcus  and  Tuckerman  (1987)  developed  another  technique 
based  on  the  spectral  methods  and  successfully  applied  it  to  the  simulation  of  a  flow 
between  concentric  rotating  spheres.  And  more  recently  Chang  and  Maxey  (1994)  have 
used  the  same  method  to  compute  the  oscillatory  flow  about  a  sphere. 


III.  GOVERNING  EQUATIONS 


A.        DERIVATION  OF  THE  EQUATIONS 

The  well  known  Navier-Stokes  equations  of  motion  of  a  viscous,  incompressible 
fluid  are; 

— +(u.V)u  =  --Vp  +  vV2u  (3.1) 

dt  p 

V.u  =  0  (3.2) 

Define  dimensionless  variables 

r    =-  u    =--  Vp    = r  V    =  — V 

a  U  v       pU2  a 

where      a:  radius  of  the  sphere 

U:  free  stream  velocity 

p  :  density  of  the  fluid 

v  :  kinematic  viscosity  of  the  fluid 

Introduce  the  dimensionless  variables  in  Eq.  (3.1) 

Uv  du       U2      .  „     .         lU2n,        U     -    . 

■+ —  (u  .V)u  =- Vp  +v-TV2u  (3.3) 


Divide  Eq.  (3.3)  by 


a2    dt        a  pa  a2 

U.v 


a2 


du       U.a     ,  _     ,        U.a„  ,       „,   „ 

-=—+ (u  \V)u   =- Vp  +vV2u  (3.4) 

dt        v  v 

Hereafter  the  superscript  *  will  be  omitted  for  nondimensional  terms.    Define  the 

U(2a) 

Reynolds  number  based  on  the  diameter,  Red  =  ,  and  take  the  curl  of  Eq.  (3.4)  to 

v 

be  able  to  eliminate  the  pressure  term  and  introduce  vorticity,  co  =  V  x  u     where  the 

vorticity  vector,  co,  is 

5 


co  =  irav+ieC0e+i<DC0kt>. 
Because  of  the  absence  of  any  azimuthal  fluid  motion  about  the  axis  of  symmetry  only  the 
O-component  of  the  vorticity  survives  and  the  transport  equation  for  the  <E>-component  of 
the  vorticity  is 


-Q),-[Vx(nx(D^  = 


2    ( 


Re 


V2- 


1 


\ 


d  v 


r2sin2e; 


co. 


(3.5) 


Convert  the  velocity  in  Eq.  (3.5)  by  defining  the  axisymmetric  Stream  function  (vj/) 


u  =  Vx 


¥ 


rsinG 


<|) ;  Azimuthal  direction 


(3.6) 


In  spherical  coordinates,  Eq.  (3.6)  is  obtained  as 


u  = 


J V 

sine  ae 


l.  - 


1      d\\f 
rsinG  3r 


where 


u  =  urir  +  ueie+u«,i0 


As  noted  earlier  u.  =  0  in  this  study 

1 


1 


(3.7) 


(3.8) 


In  the  free  stream  w  =  —  Ur2  sin2  6  or  \|/  =  —  r2  sin2  6  in  nondimensional  form. 

2  2 


By  using  the  definition  of  the  stream  function  it  can  be  shown  that 

-1    " 


C0n>  = 

2 


a2\j/   _\_  j a  (  i  d\\f 

dr2  +r2  sine  36  1  sin6  dQ 


rsin6 


(3.9) 


Define  an  operator  D    as 


_2   _a^_  sine  a  (  1 ,  d\ 
~ar2+  r2  aeUineae 


(3.10) 


to  obtain  the  azimuthal  component  of  the  vorticity; 

1     _. 


Gil>=- 


rsin6 


D> 


Call  [  V  x  (u  x  GOfcia)].^    in  Eq.  (3.5)  the  nonlinear  term  G(r,6,t) . 

By  carrying  out  the  calculations  in  spherical  coordinates  it  can  be  shown  that 


Vx(ux©0i0)  = 


1 


r3  sin  9 


a(vi/,D2\j/)    _ 


3(r,H) 


+  2D'\|/L(\|/) 


where 


u\  =  cos6 


L  = 


fi    a    i  a 


■+— 


1  —  |LL2  ar  '  r  d[i 

a(z,,z2)    az,  az2    az,  az2 

a(x,y)        dx    dy       dy    dx 


Combine  all  terms  to  get 


(3.11) 


(3.12) 


(3.13) 


a   -2  Red    1 


a(\i/,D2vj/)    _ 


d(v,\i) 


+  2D'\|fL(\|/) 


=  D> 


(3.14) 


Equation  (3.14)  is  called  the  Vorticity-Stream  function  form  of  the  Momentum 
equation  in  axisymmetric  spherical  coordinate  system.  For  computational  purposes,  which 
will  become  clear  later,  define  a  modified  stream  function  C(r,G,t),  which  is  related  to  the 
usual  Stokes  stream  function  (\|/)  by  the  relation; 

\|/  =  rCsine  (3.15) 

Also  for  any  function  f  (r,0) ,  define  a  new  operator  D2  such  that 

D2(frsine)  =  rsin9D2f 
where 


D2=V2- 


This  gives 


rz  sinz  6 

co,=-D2C  (3.16) 

u  =  Vx(CL)  (3.17) 


1         r) 

ur  = (Csin0) 

where  rsin^9  (3.18) 


The   modified   stream  function   C   is   written   as   the   sum    C  =  C  +  c,   where 
C  =  — rsinG,  (which  is  the  prescribed  uniform-flow  free  stream  condition),  and  c 

represents  the  disturbance  produced  by  the  presence  of  the  rigid  sphere. 

.  Re, 
Returning  to  Eq.  (3.5),  rearrange  by  multiplying  both  sides  with  r' 


and  the  operator  D2  is  defined  as 


2    ' 


,  Red  3co       ,  Red  ,    , 

r2-^— =  r2-^G(r,9,t)  +  r2D2co  (3.19) 

r2D2C  =  -r2©  (3.20) 

where  the  nonlinear  term 


G(r,e,t)=[Vx(iiX(D0io)].io  (3.21) 


i  d  f ,  d}       i      a 


D^=-r—  r' 


( 


+ 


r2  3rl      dr)     r2  sin2  6  39 


31  1 

sine—  -    2   .  (3.22) 

V         oQJ     r   sin   9 


B.        BOUNDARY  CONDITIONS 

The  flow  is  determined  by  numerically  solving  (3.5)  for  the  vorticity  and  (3.20)  for 
the  modified  stream  function  C  subject  to  the  boundary  conditions; 

u  =  0  at  r=  1  (3.23) 

u=U(t)    "free  stream  velocity"      as   r  -> «»  (3.24) 

Boundary  conditions  (3.23)  and  (3.24)  may  be  reformulated  in  terms  of  C  as 


dC 
C  =  0,    —  =  0  on  r=  1  (3.25) 

dr 

1 
©  =  0,  C=-U(t)rsin(0)         asr->oo  (3.26) 


The  Neumann  Boundary  Conditions  in  Equation  (3.25)  are  handled  with  the 
Green's  function  method  that  will  be  explained  later.  Homogenous  boundary  conditions 
on  the  axis  of  symmetry  (0  =  0,  n)  are  automatically  satisfied  with  the  choice  of  the  sine 
series  expansions  in  the  numerical  method  as  shown  later. 
C.        BOUNDARY  LAYER  BEHAVIOR 

The  boundary  layer  thickness  on  the  surface  of  the  sphere  can  be  calculated  with 

the  method  of  power  series  expansion  as  discussed  in,  Schlichting  (1987,pp. 235-238).  The 

body  contour  of  the  sphere  with  radius  a  is  defined  by     r(x)  =  a  sin(x  /  a)  where  x  is 

the  distance  along  the  surface  from  the  stagnation  point.  In  the  boundary-layer  limit  the 

3 
velocity  distribution  near  the  surface  of  the  sphere  is  given  by  U(x)  =  —  U^  sin(x  /  a) . 

By  introducing  the  stream-function  into  the  boundary-layer  equations  the  final 
form  of  the  transformation  is  obtained  for  axially  symmetrical  case; 


dy  d2\\f     fd\\f     ldr     \d2y     TTdU       3> 
3y  oxdy     \dx      r  dx     y  3y  dx         ay' 

with  the  boundary  conditions 

y  =  0  \|/  =  0  ,   —  =  0  (3.28a) 

dy 

y  =  oo  Jt  =  u(x)  (3.28b) 

dy 

The  potential  flow  is  defined  by  the  series  expansion  as 

U(x)  =  u1x  +  u,x3+u5x5+ (3.29) 


3U„  1U.  y     3U„a 

where    u,  = ;   u,  = r- ; T1  =  ""J 

2    a  4  a3  a  V     v 

and  the  body  contour  is  expanded  to  series  with  the  same  way; 

r(x)  =  r,x  +  r3x3  +rsx5+ (3.30) 

The  stream-function  is  represented  by  the  Blasius  series  in  analogy  with  (3.29)  & 
(3.30)  as; 

V(x,y)  =  J^-{u1xf1(TD  +  2u3x3f3(Tl)+ }  (3.31) 

By  introducing  equations  (3.29),(3.30)  and  (3.31)  into  equation  (3.27)  we  obtain  a 

set  of  differential  equations  for  the  functions  f , ,  f  3 , 

If  the  first  two  terms  of  the  expanded  series  are  concerned  the  first  equation  is 

ff=-f,f,"+i(f,,2--l)  (3.32) 

where  differentiation  with  respect  to  r|  is  denoted  by  primes.  The  boundary  conditions  are; 

r\  =  0        f,  =  f,"  =  0  (3.33a) 

Ti  =  oo        {'  =  i  (3.33b) 

10 


and  the  second  equation  is 

f,"  =  2f,'f;  -  2f,"f,  -  if.f,"  -  f,f;  -  1  (3.34) 

with  boundary  conditions 

r|  =  0         f3=f3=0  (3.35a) 

Tl  =  oo        f;=-  (3.35b) 

After  solving  the  Eq.  (3.32)  and  Eq.  (3.34)  for  f,  and  f3the  stream  function  and 

the  velocity  can  be  calculated  by  using  Eq.  (3.31).  In  the  boundary  layer  limit  the 
boundary-layer  thickness  (at  a  given  angular  location)  is  chosen  to  be  the  radial  distance 
from  the  surface  of  the  sphere  where  the  velocity  reaches  to  its  peak  value.  By  using  the 
definition  of  rj  it  can  be  shown  that 


TW=(rpe*-l)^|Red  (3.36) 

where  r    k  is  the  nondimensional  distance  where  the  velocity  reaches  its  peak  value. 

By  using  Eq.  (3.36)  the  calculated  T|peak  values  through  the  use  of  foregoing 
explanations  are  compared  with  the  corresponding  rpeak  values  which  are  obtained  from 
the  numerical  results.  The  T|peak  values  are  valid  as  Red  — >  °°  and  they  are  compared 
with  the  boundary-layer  thickness  obtained  from  the  numerical  results. 
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IV.    NUMERICAL  REPRESENTATION 


A.        COMPUTATIONAL  METHOD 

A  numerical  scheme  called  the  pseudospectral  or  collocation  method  based  on  a 
recent  paper  by  Chang  &  Maxey(1994)  which  is  in  term  based  on  a  technique  developed 
by  Marcus  &  Tuckerman  (1987)  is  used  to  solve  the  Navier-Stokes  equations.  Spectral 
methods  involve  representing  the  solution  to  a  problem  as  truncated  series  of  known 
functions.  From  this  point  of  view,  spectral  methods  may  be  viewed  as  extreme 
development  of  the  method  of  weighted  residuals  (MWR).  The  trial  and  the  test  functions 
are  the  key  elements  of  the  MWR.  This  is  also  true  for  the  spectral  methods.  The  key 
difference  lies  in  the  choice  of  the  trial  and  test  functions.  The  trial  functions  (also  called 
the  expansion  or  approximating  functions)  and  the  test  functions  (also  known  as  weight 
functions)  are  chosen  as  infinitely  differentiable  global  functions  in  spectral  methods. 
Functions  are  transformed  between  physical  space  and  spectral  space  through  the  use  of 
Fast  Fourier  Transform  (EFT)  and  Inverse  Fast  Fourier  Transform  (IFFT)  which  are 
numerically  cheap  and  quick  to  compute.  Although  they  are  not  exact  due  to  aliasing  and 
truncating  errors  they  rapidly  converges  and  few  frequency  modes  are  required  to  obtain 
accurate  results. 

A  physical  process  can  be  described  either  in  time  domain,  by  the  values  of  same  h 
as  a  function  oft,  e.g.,  h(t),  or  else  in  frequency  domain,  where  the  process  is  specified  by 
giving  its  amplitude  H,  as  a  function  of  frequency  f,  that  is  H(f),  with  -«>  <  f  <  oo.  If  h(t)  is 
real  and  even  then  H(f)  is  real  and  even  or  if  h(t)  is  real  and  odd  then  H(f)  is  imaginary  and 
odd.  The  processes  can  be  transformed  to  either  spectral  domain  or  physical  domain  by 
using  the  Fourier  transform  equations.  As  it  will  be  explained  later  in  this  chapter,  we  are 
dealing  with  odd  functions  by  using  sine  series  expansions  for  co,  c.  The  properties  of  odd 
and  even  functions  are  also  valid  for  their  derivatives  and  these  properties  will  help  to 
increase  the  computational  efficiency  in  the  calculations. 
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The  functions  which  are  introduced  into  the  governing  equations  are  presented 
both  in  spectral  space  as  a  finite  series  of  trial  functions  and  by  values  at  collocation  grid 
points  in  physical  space.  The  vorticity  CO  and  the  modified  stream  function  C  are  expanded 
as  a  Chebyshev  polynomial  in  the  radial  direction  and  as  sine  series  in  the  azimuthal 
direction.  The  calculations  in  the  radial  direction  are  done  in  physical  space  while  the 
derivatives  in  azimuthal  direction  are  obtained  in  spectral  space.  Since  each  sine  term  in 
the  expansion  satisfies  the  homogenous  0  boundary  conditions  (CO  =  0,  *¥  =  0  at  0  =  0,7t) 
exactly,  no  further  0  boundary  conditions  need  to  be  applied. 

In  general  let  a  variable,  say  co,  be  expressed  as 

Q(r,e,t)  =  5>n(z,t)SineB  (41a) 

n=l 

with 

M 

con(z,t)  =  $>mn(t)Tm(z)  (41b) 

m=0 

where  Tm(z)  is  Chebyshev  Polynomial  with  -1  <  z  <  1 

n7t 


0  = 

n      N  +  l 

n=l,2,. 

..N 

(4.2) 

z  =  Cos 

m  =  0,l,. 

..M 

(4.3) 

An  algebraic  map  is  used  to  map  the  radial  interval  1  <  r  <  r„  to  the  interval 
1  <  z  <  1  where 


r  = 


1  +  L 

V         b-zj 


lzl<l  (4.4) 


2L 

with        b  =  1  +  — (note  that  r0  =  1 ,  the  surface  of  the  sphere) 


r    —  r„ 


L  is  a  scaling  parameter  which  is  used  to  map  the  radial  interval  to  the  z-interval.  A 
finite,  but  large  outer  r.,,  was  chosen  to  avoid  regularity  problems  in  the  radial 
differentiation.    Another  mapping  parameters  is  defined  for  higher  Reynolds  numbers 
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(Red  >20).  For  small  Reynolds  numbers  a  is  taken  as  1.  Two  typical  grids  for  two 
different  parameters  are  shown  in  Figure  1 . 


Figure  1.  Typical  grids  for  a=l  a)L=8  rinf=250  M=64  N=64  b)  L=2  rinf=50  M=64  N=64 
The  sine  expansions  in  Eq.(4.1a)  and  Eq.(4.1b)  for  CO  and  C  satisfy  the 
homogenous  9  boundary  conditions  and  match  exactly  the  periodicity  and  symmetry 
conditions  in  6.  The  solution  of  the  elliptic  Eq.(3.20)  is  required  to  obtain  the  modified 
stream  function  C  from  vorticity  CO  at  any  time  level.  Following  Marcus  &  Tuckerman 
(1987)  separable  derivative  operators  Dr2  and  De2  are  defined  as 

1 


r2D2  =D2  +  • 


sin26 


D; 


(4.5) 


where 


«4 


'-*' 


D    = 


sin  9 


f 


39 


sine —  -1 

99  J 


(4.6) 


In  spectral  9-space  the  effect  of  operator  De2  at  any  fixed  radial  location  (physical  r-space) 
can  be  written  as 


D2.f(9)  = 


sin2  9 — -  +  sin9cos9 1 


39 : 


39 


I^sinCjO) 


..-If, 


j=i 


-(l  +  ij2)sin(j9)  +  ij(l  +  j)sin(29  +  j9)  +  jj(l-j)sin(29-j9) 
2  4  4 


(4.7) 
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If  we  write  this  expression  in  the  form  of  Dgf  (0)  =  Is^  If  i  we  get 


S    =  < 


(i- 2)(i-l) 


i  =  j  +  2 


-1— 


i  =  J 


(4.8) 


(i  +  2)(i  +  l) 


i  =  j-2 


0 otherwise 


This  matrix  is  NxN  if  the  representation  for  Dgf  is  truncated  to  N  Fourier  terms. 

Similarly  the  operation  of  multiplying  a  function  f  (0)  by  sin-2  0    produces  another  NxN 
matrix  which  is  called  A  matrix  where  the  only  non-zero  terms  are  ; 
'al(i)  = -i(i  + 1) i  =  j 


A;;  = 


a2(i)  =  -2i i<  j,i  +  j  =  even 


(4.9) 


B.         TIME  INTEGRATION 

Time  integration  of  the  Eq.(3.19)  and  Eq.(3.20)  are  accomplished  through  the  use 
of  an  explicit  second-order  Adams-Bashforth  scheme  for  the  nonlinear  terms  and  an 
implicit  second-order  Crank-Nicholson  scheme  for  the  viscous  and  linear  terms.  The 
calculations  are  made  in  physical  r-space  and  spectral  0-space.  If  we  denote  vorticity  co(r,t) 
as  the  vector  of  N  sine  series  then  discretized  forms  of  Eq.(3.19)  and  Eq.(3.20)  are; 


.2  Red 


co  A  -co 

t+At  t 


At 


=  r' 


Re, 


[ctGt+PGt_At]  +  [D2+A]( 


CO.  +C0 


t+At 


a 


(4.10) 


Rearrange  Eq.(4.10)  and  get ; 


D2+A-r" 


Re, 

At 


©,+*  = 


D2+A  +  r2R6d 


At 


cot-r2Red[cxGt-pGt_At]  (4.11) 


3  1  r  -, 

where    a  =  — ;p  =  -—     and    G(r,0,t)  =  [Vx(uxco)Ji(> 

Equation  (3.20)  is  also  discretized  by  using  the  same  method  as  ; 
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[Dr2+A]c(r,t  +  At)  =  -r2co(r,t  +  At)  (4.12) 

The  radial  derivatives  are  evaluated  by  collocation  methods  and  expressed  as 
matrix  operations  on  the  vector  of  function  values  at  the  grid  points  in  physical  r-space. 
For  (M+l)  collocation  points,  we  introduce  the  (M+l)x  (M+l)  matrix  operation,  that  is 

Re, 


Dw'(n)  =  D,  + 


a(D- 


At 


'dr2 


I 


(4.13) 


where,  I  is  the  identity  matrix.  The  operator  [Dr  -r  ( Red/At)+A]  in  Eq.(4.1 1)  may  now  be 
written  in  block  matrix  form  as 


D(1)(l) 

0 

a2(l)I 

0 

0 

D(1)(2) 

0 

a2(2)I     ... 

0 

0 

D(1)(3) 

0 

0 

0 

0 

D(,)(4)     ... 

•  •  • 

•  •  • 

•  •  • 

... 

(4.14) 


This  matrix  acts  on  the  vector  of  a  length  of  (M+l)  x  N  for  vorticity  coefficients. 
Eq.(4.11)  and  Eq.(4.12)  are  upper  triangular,  block  matrix  problems  in  physical  r-space 
and  spectral  9-space  that  can  be  solved  by  any  solver  with  the  Dirichlet  boundary 
conditions  discussed  below.  The  Eq.(3.32)  is  treated  in  the  same  manner  and  leads  to  an 
upper  triangular,  block  matrix  problem  with  each  diagonal  term  on  row  n  replaced  by 
Dr2+al(n)I. 
C.        NONLINEAR  TERMS 

If  one  opens  the  nonlinear  term  G(r,0,t)  in  spherical  coordinates  gets  six  terms; 


1  dCQ0  dC 
r    3r    90 


cot(9)     da 


3r 


1  to,  dC 

r   36    3r 


r2  ae 


cot(8)     ac 


(4.15) 


C  9<P+ 

r2  ae 


The  nonlinear  terms   in  Eq.(4.15)   are  handled  by  transforming   the   spectral 
coefficients  to  physical  space  first  and  transforming  the  result  of  the  products  to  spectral 
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space  after  the  multiplication.  Basically  the  derivatives  are  obtained  in  spectral  space  and 

the  products  are  done  in  physical  space  and  finally  the  result  is  transformed  to  spectral 

space. 

D.        GREEN'S  FUNCTION  METHOD 

In  order  to  enforce  the  radial  boundary  conditions  in  Eq.(4.11)  and  Eq.(4.12), 
Green's  function  method  has  been  developed  which  gives  the  correct  boundary  conditions 
in  Eq.(3.25)  and  Eq.(3.26)  and  allows  the  surface  vorticity  to  develop  naturally.  If  co  and 
c  consist  of  two  parts,  homogenous  and  particular  solutions,  we  write 

N  N 

(D  =  COp  +X^jSj  and  c  =  cp+X?1jCj         J  =  1,2, N  (4.16) 

j=i  j=i 

If  we  introduce  Eq.(4.16)  into  Eq.(4.11)  and  Eq.(4.12),  we  find  the  homogenous  parts 
with  corresponding  boundary  conditions  as  follows; 

E2coj(r,6)  =  0  (4.17)     H'c^-r'cOj  (4.18) 

Sj(r  =  l,e.)  =  5ij  ffJ(r  =  l,B1)  =  0 

®.(r  =  oo,e.)  =  0  cj(r  =  oo,e.)  =  0 

The  particular  parts  with  corresponding  boundary  conditions  are  ; 

E2cop(r,6,t  +  At)  =  R(r,e,t)  (4.19)      H2cp(r,0,t  +  At)  =  -r2cop(r,6,t)       (4.20) 

C0p(r  =  l,9)  =  0  Cp(r  =  l,e)  =  -Clr=1 

©  (r  =  oo,6)  =  0  c  (r  =  oo,e)  =  0 

p 

E2  and  H2  are  obtained  from  block  matrix  defined  as  in  Eq.(4. 14) 

9c  dC 

To  find  A,  we  enforce  -r-lr_,  =  — -r— L,  (4.21) 

dr  dr 

If  we  substitute  Eq.(4.16)  into  Eq.(4.21),  we  get 


3c  n  3c  3c 

+  I^7    I    *,~f    I  (4.22) 


dr  r=i,e.      j=i  3r  r=1,6|  3r  r=1,e, 
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Let 


3c      3c, 


3r      3r 


Jr=i,ei 


3c. 
=  Hi        and       -J-    |    =A, 

OT    r=l,e, 


(4.23) 


finally  from  Eq.(4.21),Eq.(4.22)  and  Eq.(4.23),  we  find     Hj  =XAij^j  and  solve  for 


j=i 


X-    asX  =  A"1!!.  These  X  values  can  now  be  used  in  Eq.(4. 16)  to  arrive  at  the  values  of  co 
andc. 
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V.  DISCUSSION  OF  THE  RESULTS 

The  results  are  going  to  be  discussed  in  three  parts.  The  drag  coefficients 
calculated  by  following  the  explanation  below  for  various  Reynolds  numbers  are  compared 
with  the  results  of  Chang  &  Maxey  (1994)  and  the  streamlines  contours  are  compared 
with  the  experimental  flow  visualization  pictures  and  finally  the  boundary-layer  behavior  in 
the  boundary-layer  limit  is  discussed  by  using  the  numerical  results. 
A.        EVALUATION  OF  THE  DRAG 

By  integrating  the  contributions  from  the  stress  tensor  C-  over  the  sphere  surface 

the  resultant  fluid  force  on  the  sphere  exerted  by  the  moving  fluid  is  found  as; 

Fj  =|<y^nJdS  (5.1) 

where  n  is  the  unit  outward  normal.  There  is  no  lift  force  because  the  flow  is  axisymmetric 
and  the  force  F  acts  parallel  to  the  free-stream  velocity  U.  The  components  of  the  force 
can  be  expressed  in  terms  of  the  spherical  polar  coordinates  as  ; 

F,  =  j(a„  cos(9)-are  sin(6))dS  (5.2) 

The  normal  component  of  the  stress  a^  is 

crT=-p  +  2pv-^L  (5.3) 

The  no  slip  boundary  conditions  in  Eq.(3.23)  and  the  condition  of  incompressible  flow  in 
Eq.(3.2)  causes  the  dujdr  to  vanish  on  the  surface  of  the  sphere  and  the  only 
contribution  from  Gn  is  through  pressure.  The  shear  stress  ar8  is 


<*r9  =  2PV 


r  9  fue^     1  3ur 


2  dr I   r  J     2r  39 
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(5.4) 


Equation  (5.3)  and  Eq.(5.4)  are  in  dimensional  form.  By  using  the  non-dimensional 
variables  the  drag  force  Cd  is  defined  as 

Cd=F,/7tpa2U^  (5.5) 

The  drag  force      Cdmay  be  split  into  two  separate  parts,  Cffor  the  frictional 

component  due  to  shear  stress  and      Cpfor  the  pressure  component.  The  frictional 

component  is  calculated  from  Eq.(5.2)  and  Eq.(5.4)  as 

4    f 
Cf  =-  —  lco(r=l,0)sin2(0)d0  (5.6) 

Re  0 

The  pressure  component  Cp  in  non-dimensional  form 

K 

Cp  =  -Jp(r  =  l,0)sin(20)d0  (5.7) 

o 

The  pressure  variation  over  the  surface  of  the  sphere  can  be  determined  once  the 

vorticity  field  is  known  : 

p(r  =  l,e)  =  p(l,7c)-  — J  — (rco)  I  d0'  (5.8) 

Re  J0  dr         r=, 

A   non-dimensional   pressure   coefficient    k(6)    can   be    defined    to    make    the 
calculations  easier; 

k(9)  =  2p(r  =  l,0)  (5.9) 

The  pressure  at  the  stagnation  point,  G  =  n ,  can  be  evaluated  from  the  radial 
momentum  equation  as; 

k(7c)  =  2p(r  =  l,7C)  =  l  +  —  ]-—  |  dx  (5.10) 

Re  -J  r  30  e=n 

B.       STEADY  FLOWS 

In  order  to  verify  the  numerical  method  and  to  provide  a  comparison,  simulations 
for  steady,  unidirectional  flow  were  carried  out  for  Re  <  1 04 .  The  parameters  which 
were  used  in  the  calculations  are  listed  in  Table  1 . 
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For  all  calculations  the  flow  is  initially  at  rest  and  then  a  uniform,  constant  external 
flow  is  introduced  in  the  direction  of  0  =  n  and  the  simulation  continued  till  the  solution 
has  converged  to  steady  state.  Two  different  grids  were  used  in  the  calculations.  A  spatial 
discretization  of  64x64  grid  points  was  used  for  Re  <  10  and  a  relatively  finer  mesh  of 
128x128  grid  points  was  used  for  Re>  20.  For  Re  >  20  a  mapping  parameter  "a"  is 
defined  for  finer  resolution  near  the  surface.  The  vorticity  diffuses  over  a  large  region  at 
low  Reynolds  numbers  compared  to  the  sphere  radius  while  it  is  convected  downstream 
with  a  more  complex  structure,  which  needs  increased  resolution. 


Re 

L 

a 

M 

N 

At 

r„ 

0.1 

8 

64 

64 

0.05 

250 

0.2 

8 

64 

64 

0.05 

250 

0.3 

8 

64 

64 

0.05 

250 

0.4 

8 

64 

64 

0.05 

250 

0.5 

8 

64 

64 

0.05 

250 

0.6 

8 

64 

64 

0.05 

250 

0.7 

8 

64 

64 

0.05 

250 

0.8 

8 

64 

64 

0.05 

250 

0.9 

8 

64 

64 

0.05 

250 

1.0 

8 

64 

64 

0.05 

250 

5 

2 

64 

64 

0.02 

50 

10 

2 

64 

64 

0.02 

50 

20 

1.5 

1.5 

128 

128 

0.02 

50 

30 

1.5 

1.5 

128 

128 

0.01 

40 

40 

1.5 

1.5 

128 

128 

0.01 

30 

56.5 

1.5 

1.5 

128 

128 

0.01 

30 

104 

1.5 

1.5 

128 

128 

0.005 

25 

Table  1 .  Parameters  used  in  the  calculations 
The  present  results  have  been  summarized  in  Table  2  with  the  previous  results 
which  were  calculated  by  Chang  &  Maxey  (1994).  As  it  can  be  seen  in  Table  2  there  is  a 
general  agreement  between  the  present  results  and  the  previous  results.  The  separation 
angle,  0S,  and  the  ratio  of  bubble  length  (s)  to  the  diameter  of  sphere  (d)  are  plotted  later 
in  this  chapter  to  compare  them  with  the  results  which  were  obtained  experimentally  by 
Taneda(1956). 
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Figure  2  shows  the  variation  of  vorticity  and  streamlines  in  the  vicinity  of  the 
sphere  at  Re  =  0. 1 .  As  it  can  be  noticed  the  vorticity  is  diffusing  over  a  large  region  and 
very  slight  asymmetry  about  0  =  n  1 2  is  observable. 
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Figure  2.  (a)  Lines  of  constant  vorticity  for  Re  =  0.1  Aco  =  0.1  (b)  Streamlines  A\|/  =  0.3 

Figure  3  shows  the  constant  vorticity  and  streamlines  contours  for  Re  =  5  .The 

asymmetric  region  is  getting  more  observable  about   6  =  %  1 2      as  Reynolds  number 

increases. 
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Figure  3.  (a)  Lines  of  constant  vorticity  for  Re  =  5  Aco  =  0.1 5  (b)  Streamlines  A\j/  =  0.25 

In  Figure  4  the  constant  vorticity  and  streamlines  contours  are  plotted  for  Re  =  40. 

A  separation  region  is  clearly  observable  on  both  plots.  The  recirculation  in  the  separation 

region  is  also  noticeable  with  the  asymmetry  over  the  sphere. 
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Figure  4.  (a)  Lines  of  constant  vorticity  for  Re  =  40  Ago  =  0.1  (b)  Streamlines  A\j/  =  0.05 
for  ' — ',  A\j/=  0.005  for  '-.-.-.-',  Ay  =  0.00004  for  the  separation  region. 

Figure  5  magnifies  the  separation  region  over  the  constant  streamline  contours  for 
Re  =  40.  There  is  a  distortion  of  the  vorticity  field  caused  by  advection  in  the  flow  and  a 
separation  is  clearly  visible  with  a  small  region  of  negative  streamlines  in  recirculation. 
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Figure  5.  Constant  streamline  contours  for  Re  =  40.  A\}l  =  0.00004  in  the  separation 
region  ,  A\|/=  0.005  for  '-.-.-.-',  A\j/=  0.2  for  ' ' 
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Figure  6  shows  the  constant  vorticity  and  streamlines  for  Re  =  56.5.  The 
separation  and  the  recirculation  inside  the  separation  region  is  clearly  visible.  There  is  a 
noticeable  change  in  the  wake  length  as  the  Reynolds  number  is  increased  from  Re  =  40  to 
Re  =  56.5. 


0  2  -2  0  2 

Figure  6.(a)  Constant  vorticity  contours  for  Re  =  56.5  Aft)  =  0.3(b)  Streamline  contours 
Ay  =0.005  for  * — ',  Ay  =  0.15  for  '-.-.-.-',  Ay  =  0.0003  for  the  separation  region. 

The  constant  vorticity  and  streamline  contours  for  Re  =  104  are  plotted  in  Figure 
7.  The  separation  region  and  recirculation  inside  the  separated  flow  is  clearly  visible.  The 
increase  in  the  wake  length  is  also  noticeable  for  this  particular  Reynolds  number.  It  is  also 
noticeable  from  the  foregoing  plots  that  the  vorticity  diffuses  to  smaller  region  as 
Reynolds  number  increases. 
4p ! ! 1  ^t  =  :  =  ^--\  = 
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Figure  7.  (a)  Lines  of  constant  vorticity  for  Re  =  104  Aco  =  0.4  (b)  Streamlines  Ay  =0.005 
for  ' — ',  Ay  =  0.15  for  '-.-.-.-',  Ay  =  0.0009  for  the  separation  region  for  Re  =  104. 


27 


Figure  8  magnifies  the  streamlines  for  Re  =  56.5.  The  flow  visualization  picture 
for  same  Reynolds  number  is  also  presented  in  Figure  9.  It  is  clearly  observable  that  there 
is  a  similarity  between  Figure  8  and  Figure  9.  The  wake  length,  separation  location  and  the 
nature  of  the  streamlines  in  Figure  8  and  Figure  9  look  alike  in  both  figures  for  the 
numerical  results  and  flow  visualization  picture. 

Figure  10  magnifies  the  streamlines  plot  for  Re  =  104  to  present  the  similarity  of  it 
with  the  flow  visualization  picture  in  Figure  11  for  the  same  Reynolds  number.  The 
extension  of  the  recirculating  wake  to  a  full  diameter  in  the  downstream  is  clearly  visible  in 
Figure  10  and  Figure  11.  The  flow  is  still  laminar  at  this  Reynolds  number. 

The  calculated  separation  angles  and  the  ratio  of  the  wake  lengths  to  the  diameter 
of  the  sphere  are  presented  in  Figure  12  respectively  along  with  the  experimental  data 
which  were  obtained  by  Taneda  (1956).  There  is  a  good  agreement  between  the  present 
results  and  Taneda' s  data.  Taneda  has  also  extrapolated  his  data  and  conjectured  the 
critical  Reynolds  number  of  the  "onset  of  separation"  as  Rec  =  24  .On  the  other  hand,  the 

earlier  work  of  Nisi&Porter  (1926), which  was  specially  designed  to  detect  the  onset  of 
separation  with  ultramicroscopic  methods,  reported  separation  at  Re  =  10  and  indicated 
that  the  onset  of  separation  in  an  infinite  fluid  is  about  Rec  =  8 . 

An  estimate  of  Re  =  20  is  given  by  the  present  results  as  the  Reynolds  number  for 
the  onset  of  separation.  This  is  good  agreement  with  the  value  of  Rec  =  20.7  given  by 

Chang&Maxey  (1994),  the  value  of  Rec  =  20.5  given  by  Dennis&Walker  (1971)  and  the 
value  of  Rec  =  20given  by  both  Le  Clair  et  al.(1970)  and  Lin&Lee  (1973).  Others  have 
found  separation  to  first  occur  at  different  Reynolds  numbers,  for  example  Rimon&Cheng 
(1969)  found  separation  as  low  as  Rec  =10. 

Taneda' s  data  also  suggests  that  the  length  of  the  separation  bubble  varies  linearly 
with  the  logarithm  of  the  Reynolds  number.  The  present  results  are  consistent  with  such  a 
result  and  they  remark  that  over  this  limited  range  a  simple  linear  relationship  with 
Reynolds  number  is  a  good  approximation. 
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Figure  8.  Streamlines  A\j/  =0.005  for  ' — ',  A\|/=  0.15  for  '- 

separation  region  for  Re  =  56.5 
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Figure  9.  The  flow  visualization  picture  of  streamlines  for  Re  =  56.5  (  Archives  de 
l'Academie  des  Sciences  de  Paris.  Payard  &  Coutanceau  1974  ) 
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Figure  10.  Streamlines  A\\f  =0.005  for  ' — \  A\j/=  0.15  for  '- 

separation  region  for  Re  =  104. 
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Figure  11.  Visualization  picture  for  Re  =  104  by  a  thin  coating  of  condensed  milk  on  the 
sphere,  which  gradually  melts  and  is  carried  into  the  stream  of  water.  Taneda  1956. 
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Figure  12.(a)  Separation  angle  0S  vs.  Log  Re  (b)  The  ratio  of  bubble  length  to  diameter 

(s/d)  vs.  Log  Re  '— '  Taneda  (1956) 

C.       BOUNDARY  LAYER  BEHAVIOR 

Once  the  stream  function  was  known  for  the  particular  flow  the  velocity  over  the 
sphere  can  be  calculated  by  using  the  Eq.(3.7).  The  velocity  in  0  direction  ue  which  is 
parallel  to  the  surface  of  the  sphere  can  be  calculated  by  using  the  relation 

1      3\|/ 

u°= — nnr  (5-n) 

rsinQ  or 
We  also  know  from  potential  flow  theory  that  the  stream  function  is  given  by 

\|/  =  -(r2-i)sin2e  (5.12) 

By  using  the  relation  in  Eq.(5.1 1)  uecan  be  calculated  from  Eq.(5.12)  for  potential  flow 
as 

T2r3  +  0 


U9   = 


2r3 


sinG 


(5.13) 


Figure  13  shows  the  ue/sin0  values  which  are  obtained  from  Eq.(5.13)  for 
potential  flow  along  with  the  numerical  results  for  various  Reynolds  numbers  from  this 
study.  It  is  clearly  observable  that  the  velocity  profile  matches  with  the  potential  flow  at 
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higher  Reynolds  numbers.  The  velocity  profile  for  a  small  Reynolds  number  is  also 
presented  in  Figure  13  for  comparison,  and  shows  that  the  velocity  reaches  its  free  stream 
value  over  a  much  larger  distance  in  this  Stokes  flow  regime. 


i-5rr 


Figure  13.  The  ue  /sinG  profiles  for  varying  Reynolds  numbers  (0  =  n  I  4 from  the  nose) 
' '  represents  potential  flow,  solid  lines  represent  numerical  results. 

The  approach  to  boundary-layer  behavior  is  now  examined  by  comparing  the 

boundary-layer  thickness  from  numerical  results  with  values  of  T\  peak  from  boundary-layer 

analysis  as  discussed  earlier  in  Chapter  III.C.  Figure   14  shows  the   T)    k    values  at 

corresponding  Reynolds  numbers  along  with  the    T|peak  value  as   Re— »«>  for  a  particular 

angle. 

The  T|    .  values  for  corresponding  Reynolds  numbers  can  be  modeled  by  assuming 
a  model  such  as 

r|  =  C.+ — -  (5.14) 

1     Rem 

where  C ,  is  the  Tj  peak  value  as  Re— »©o  ,  C  2  and  m  are  the  constants. 
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Figure  14.  T|peak  values  vs.  Log  Re  ' — '  represent  the  T|peak  value  as  Re— >°° 

By  using  the  T|peak  values  at  Re=30,  Re  =  56.5  and  Re  =  104  ,  C2  and  m  can  be 
calculated  from  Eq.(5.14)  as  C2  =  2.464  and  m  =  0.182  to  obtain  the  model  as 

2.464 


T1  =  C,+ 


Re 


0.182 


(5.15) 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  numerical  results  have  showed  a  good  agreement  with  the  other  results  which 
were  obtained  by  using  different  numerical  methods.  The  pressure  coefficient,  C  ,  which 

is  one  of  two  components  of  the  total  drag  coefficient,  Cd ,  couldn't  be  resolved  for  higher 
Reynolds  numbers  through  the  use  of  present  method.  But  the  calculated  skin  friction 
coefficients,  Cf ,  and  the  stagnation-point  pressure  k(n)  have  been  resolved  successfully. 

A  fine  mapping  for  collocation  points  is  very  important  to  obtain  more  accurate 
results.  The  mapping  parameters  r^ ,  OC  and  L  should  be  chosen  in  a  way  to  provide 

enough  collocation  points  in  the  boundary  layer  and  recirculation  regions.  A  finer  mapping 
can  increase  the  accuracy  by  allowing  finer  resolution  in  these  large  gradient  regions. 

The  Green's  function  method  that  provides  the  natural  development  of  vorticity 
along  with  the  Neumann  boundary  conditions  should  not  be  considered  as  a  trivial  detail 
and  must  be  used  to  have  a  stable  convergence  process. 

A  different  FFT  technique  can  be  used  which  does  not  restrict  the  number  of  grid 
points  to  a  power  of  two.  This  could  reduce  the  memory  requirement  by  allowing  a 
relaxation  in  the  choice  of  grid  points. 

Very  good  agreement  has  been  obtained  with  the  previous  calculations  and 
available  experimental  flow  visualization  pictures.  This  study  is  at  a  stage  where  it  can 
now  be  extended  to  explore  higher  Reynolds  numbers  with  suitable  turbulent  models. 
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APPENDIX  A.  PROGRAM  STRUCTURE 


In  the  first  part  of  the  program,  problem  parameters  were  defined  such  that  Re,  dt, 
L,  r_inf,  M,  N.  Calculations  which  don't  depend  on  time  were  also  performed  in  this  part. 
In  order  to  gain  some  additional  memory  Dl,  D2  and  D3  matrices  in  (4.1 1)&(4.12)  were 
calculated  with  the  function  D_Bmat  and  saved  in  terms  of  diagonal  elements  only.  Non- 
diagonal  elements  of  these  matrices  were  included  separately  in  each  iteration. 

Calculations  showed  that  the  matrix  Ml  (derivative  operator  in  r  direction)  and 
consequently  the  matrix  M3  and  the  block  matrices  which  are  formed  by  the  diagonals  of 
Dl,  D2  and  D3  are  ill-conditioned.  To  avoid  the  numerical  errors  in  solving  the  systems  of 
equations,  Singular  Value  Decomposition  (S VD)  method  was  used.  SVD  of  the  matrices 
and  calculations  of  w  and  c  that  are  used  in  Green  Function  method  were  performed  in  this 
part.  The  built-in  function  SVD  in  MATLAB,  was  used  for  decompositions  of  the 
matrices  in  the  systems  of  equations. 

The  calculation  part  basically  consists  of  four  steps.  In  first  step,  A  ^matrix  was 

calculated  from  eqn.  (4.22)  out  of  the  time  loop.  In  second  step,  inside  the  time  loop  the 
particular  solutions  of  w  and  c  were  found  from  eq.  (4.17)  and  (4.18).  In  third  step,  H  that 
is  used  in  Green's  Function  method,  was  calculated  and  finally  homogenous  and  particular 
solutions  are  combined  to  get  w  and  c. 

Most  of  the  calculations  were  performed  in  spectral  domain.  Only  nonlinear 
operations  were  done  in  physical  domain.  Nonlinear  terms  were  calculated  with  the 
functions  nonlin  for  w-c  equations  . 

A  convergence  criterion  isn't  used  in  the  program  to  avoid  the  wrong  results  which 
are  caused  by  the  over-shooting  for  high  Reynolds  numbers.  Instead,  the  results  are 
checked  during  the  calculations  and  program  was  terminated  after  seeing  the  steady  state. 

Solsvd  :  This  function  is  used  to  solve  equations  after  SVD  in  Matlab.  A  limit, 
which  is  defined  as  a  number  smaller  than  the  values  of  S  matrix  in  SVD,  is  used  in  this 
function  to  avoid  the  numerical  errors  due  to  those  ill-conditioned  matrices.  This  is 
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avoided  by  forcing  those  singular  values  which  are  smaller  than  the  defined  limit,  to 
zero.  10"14  is  used  in  all  calculations. 

Solsvdl  :  Since  this  function  was  used  for  solving  homogenous  parts  in  Green's 
Function  method  and  performed  once  before  the  time-loop,  no  limit  was  used  in  solsvdl. 

Phy  :  This  function  was  used  for  transformations  of  the  variables  from  spectral 
domain  to  physical  domain 

Operat  :  Derivatives  in  0  direction  in  spectral  domain  was  performed  with  the 
function,  operat.  Since  these  terms  were  used  only  in  nonlinear  terms  the  output  from  this 
function  is  in  physical  domain. 

L       :  Mapping  parameter 

Re      :  Reynold's  Number 

dt     :  Time  interval  between  iterations 

r_inf :  Outer  boundary 

step  :  Number  of  iterations 

N       :  Number  of  grid  points  in  0  direction. 

M       :  Number  of  grid  points  in  r  direction. 

xcor,  ycor  :  Cartesian  coordinates  of  grid  points  with  respect  to  center  of  sphere, 
(used  for  presentation  of  the  results) 

Ml      :  Derivative  operator  in  r  direction. 

Dl,  D2,  D3  :  Matrixes  defined  as  in  (4.1 1)&(4.12) 

oo        :  Vorticity. 

c         :  Potential  function  as  in  (3.15a)&(3.26) 

wtil,  ctil :  Homogenous  part  of  the  co  and  c  defined  as  in  Green  Function. 

Aij,  H,  landa  :  Variables  defined  as  in  Green's  function  method  in  chapter  IV. 

wp,  cp  :  Particular  solution  of  co  and  c  defined  as  in  Green  Function. 

Cf       :  The  frictional  component  of  drag 

kpi        :  Non-dimensional  pressure  coefficient  at  0  =  n 
Cp        :  Pressure  component  of  drag 
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APPENDIX  B.  PROGRAM  CODES 


%%  define  program  parameters 

%  Table  1  list  these  parameters  for  every  Reynolds  numbers 

clear  all 

global  E  Ml  M2M3 

L=8;Re=l;dt=0.02; 

r_inf=250;step=  10000; 

%%  define  grid  points 

N=32;n=[0:N-l]'; 
teta=n.*pi./N; 

M=32;m=[0:M]'; 

z=cos(pi.*m./M); 
b=l+2*L/(r_inf-l);  %  map=l 
r=l+L.*(l+z)./(b-z);%  map=l 
global  NM 
xcor=cos(teta)*r'; 
ycor=sin(teta)*r'; 

E=emat(M); 

disp('E  matrix  is  done') 

%%  define  Ml, M2,M3 

dz_dr=(b.*(L+r-l)-b.*(r-l)+L)./(L+r-l).A2;  %  map=l 

B=diag(dz_dr); 

M1=B*E; 

R=diag(r);global  R  Re  dt 

M2=R*M1; 

M3=R*(2.*eye(M+ 1  )+M2)*M  1 ; 

%%  define  Dl  D2  D3 

[Dl,D2,D3]=D_Bmat(M,N); 

%  define  initial  w  c  wp  cp  u 

duml=(-1.5./r.A2)*sin(teta'))'; 

dum2=[duml;zeros(l,M+l);-l.*flipud(duml(2:N,:))]; 

dum3=fft(dum2); 
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w=imag(dum3(l:N,  :))];w=w'; 
duml=(-0.75+0.25./r.A2)*sin(teta'))'; 
dum2=[duml;zeros(l,M+l);-l.*flipud(duml(2:N,:))]; 
dum3=fft(dum2);c=imag(dum3(  1  :N,:  ))];c=c' ; 
wp=zeros(M+l  ,N);cp=wp; 

%  define  kmat  (for  teta  der.)  and  cotteta  (for  nonlinear  terms) 

dum=ones(  1  ,M+ 1 ) ; 
kmat=n*dum; 

cotteta=[0;cot(teta(2:N))]  *dum; 
global  kmat  cotteta  teta 

%  svd  decomposition  for  Dl  D3  and  define  rsquare  (nonlinear) 

Dlu=zeros(N*(M+l),M+l); 

Dll=Dlu;Dlp=Dlu;D31=Dlu;D3u=Dlu;D3p=Dlu; 

rsquare=zeros(M+l  ,N); 

fori=l:N 
bas=(i-l)*(M+l)+l; 
son=bas+M; 
rsquare(:,i)=r.A2; 

Dl(bas,:)=[l  zeros(l,M)];  %  modify  Dl  for  w  @  r=inf 
Dl(son,:)=[zeros(l,M)  1];  %  modify  Dl  for  w  @  r=l 
[Dll(bas:son,:),Dlu(bas:son,:),Dlp(bas:son,:)]=svd(Dl(bas:son,:)); 
D3(bas,:)=[l  zeros(l,M)];  %  modify  D3  for  c  @  r=inf 
D3(son,:)=[zeros(l,M)  1];  %  modify  D3  for  c  @  r=l 
[D31(bas:son,:),D3u(bas:son,:),D3p(bas:son,:)]=svd(D3(bas:son,:)); 

end 

global  rsquare 

%  define  c_  c_teta  c_r  u_  u_teta  u_r 

c_=(0.5.*sin(teta))*r'; 
c_teta=0.5.*cos(teta))*r'; 
c_r=-(0.5*sin(teta))*ones(  1  ,M+ 1 ); 
global  c_  c_teta  c_r 

%  define  c  boundary  at  r=l 

dum4=-0.5.*sin(teta); 

dum5=[dum4;0;-l  .*  flipud(dum4(2:N,:))]; 

dum6=fft(dum5); 

cbr  1 =imag(dum6(  1  ,N))' ; 
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%  begin  green  function  wtil  ctil 

wtil=zeros(M+l  ,NA2); 
ctil=wtil; 

forjj=l:N 
RR=zeros(M+l,N); 
duml=[zeros(l,jj-l)  1  zeros(l,N-jj)]; 
dum2=[duml  0  -l.*fliplr(duml(2:N))]; 
dum3=fft(dum2); 

dum4=imag(dum3(  1  :N));if  jj==  1  ;dum4=real(dum3(  1  :N));end; 
RR(M+l,:)=dum4; 
wtil(:jj*N)=solsvdl(Dll((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),. 

Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

RR(:,N)); 

fori=N-l:-l:l 
sum=zeros(M+ 1,1); 
ii=(jj-l)*N+i; 
forj=i+l:N 

if  rem(i+j,2)=2 
sum=sum+(-2*(i- 1  )).*wtil(:,(jj- 1  )*N+j); 

end 
sum(l)=0; 
sum(M+l)=0; 
end 
RR_=RR(:,i)-sum; 

wtil(:,ii)=solsvdl(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
RR_); 

end 

tempwtil=(-l.*rsquare).*wtil(:,(jj-l)*N+l:jj*N); 
temp  wtil(M+ 1 , :  )=zeros(  1  ,N) ; 
tempwtil(  1  ,:)=zeros(  1  ,N) ; 
ctil(:,jj*N)=solsvdl(D31((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),.. 

D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

tempwtil(:,N)); 
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fori=N-l:-l:l 
ii=(jj-l)*N+i; 
sum=zeros(M+ 1,1); 
forj=i+l:N 

ifrem(i+j,2)=2 
sum=sum+(-2*(i-l)).*ctil(:,(jj-l)*N+j); 

end 
sum(l)=0; 
sum(M+l)=0; 
end 
wtil_=tempwtil(:,i)-sum; 

ctil(:,ii)=solsvdl(D31((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),. 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
wtil_); 
end 
end 

%  find  Aij  matrix 

forjj=l:N 

blok=ctil(:,(jj- 1  )*N+1  :jj*N); 

duml=phy(blok')'; 

dum2=(Ml*duml)'; 

dum3=[dum2;  zeros(l,M+l) ;  -l.*flipud(dum2(2:N,:))]; 

dum4=fft(dum3);dum5=imag(dum4(  1  :N,:));dum5=dum5'; 

ctilr_l(jj,  1  :N)=dum5(M+l ,:); 
end 

Aij=ctilr_l'; 

%  find  svd  of  Aij  for  landa  solving 

Aij(l,:)=[l  zeros(l,N-l)]; 
[Aiju,Aijs,Aijv]=svd(Aij); 

%  begin  time  loop 

for  ops=l:  step 

if  ops==l  ;prc=c;prw=w;pru=u;end 
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%  findRl 

Rl=zeros(M+l,N); 

Rl(:,N)=D2((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:)*w(:,N); 
fori=N-l:-l:l 
summ=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)=2 

summ=summ+(-2*(i-l)).*w(:,j); 
end 
end 

Rl(:,i)=D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)*w(:,i)+summ; 
end 
R1=-1.*R1; 

%    findR2 

R2a=nonlin(w,c); 
R2b=nonlin(pr  w  ,prc) ; 
R2=3/2.*R2a-l/2.*R2b; 
R2=(- 1  *Re).*rsquare.*R2; 
Rtot=Rl+R2; 
prw=w;prc=c; 

%%  green  function  solution 

%  find  new  wp 

%%  apply  BC  for  wp 

Rtot(  1  ,:)=zeros(  1  ,N);    %  BC  for  wp=0  at  r_inf 
Rtot(M+l,:)=zeros(l,N);  %  BC  for  wp=0  at  r_l 

wp(:,N)=solsvd(Dll((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),.-. 
Rtot(:,N)); 

fori=N-l:-l:l 
sum=zeros(M+ 1,1); 
forj=i+l:N 

if  rem(i+j,2)=2 

sum=sum+(-2*(i-l)).*wp(:,j); 
end 
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sum(l)=0; 
sum(M+l)=0; 
end 
Rtot_=Rtot(:,i)-sum; 

wp(:,i)=solsvd(Dll((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)f... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
RtotJ; 

end 

%  solve  for  new  cp 

tempwp=- 1  .*rsquare.*  wp; 

tempwp(M+l,:)=cbrl;  %  BC  for  cp=-c_r  at  r_l  zero  for  pure  spin 

tempwp(l,:)=zeros(l,N);  %  BC  for  cp=0  at  r_inf 

cp(:,N)=solsvd(D31((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
D3p((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
tempwp(:,N)); 

fori=N-l:-l:l 
sum=zeros(M+l,l); 
forj=i+l:N 

ifrem(i+j,2)=2 
sum=sum+(-2*(i-l)).*cp(:,j); 

end 
sum(l)=0; 
sum(M+l)=0; 
end 
w_=tempwp(:,i)-sum; 

cp(:,i)=solsvd(D31((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
w_); 
end 

%  find  landa 
%  find  H 

duml=phy(cp')'; 
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dum2=(Ml*duml)';  %  NxM+1 

dum3=- 1  .*c_r-dum2; 

dum4=[dum3;zeros(l,M+l);-l.*flipud(dum3(2:N,:))]; 

dum5=fft(dum4); 

dum6=imag(dum5(  1  :N,:)); 

dum6=dum6';  %  M+lxN 

H=dum6(M+l,:);H=H';  %  Nxl 

H(1)=0; 

landa=solsvd(Aiju,Aijs,Aijv,H); 

%  find  complete  w  &  c 

dum  1  c=zeros(M+ 1  ,N)  ;dum  1  w=dum  1  c; 
forjj=l:N 

blokc=ctil(:,(jj- 1  )*N+ 1  :jj  *N); 

blokw=wtil(:,(jj-l)*N+l:jj*N); 

dum2c=landa(jj ) .  *blokc ; 

dum2  w=landa(jj ) .  *blokw ; 

dum  1  c=dum  1  c+dum2c; 

dum  1  w=dum  1  w+dum2  w ; 
end 

w=wp+dumlw; 
c=cp+dumlc; 

res=phy(w'); 
resp=phy(prw'); 
resc=phy(c'); 
resc=resc+c_; 

%  Calculate  Cf  for  every  iteration 
duml=res(:,M+l); 
dum2=(sin(teta)).A2; 
dum3=dum  1  .*dum2; 
cf(ops)=-4/Re*trapz(teta,dum3); 

%  Calculate  k(7t)  for  every  iteration 

duml=w'; 
fori=l:M+l 

sum(i)=0; 

forj=l:N 

f=(a-l)*duml(i,j).*(-l.Aa-D))-/N; 
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sum(i)=sum(i)+f; 
end 
end 

dum2=(  1  ./r).*sum; 
kpi(ops)=(l+(8/Re)*trapz(r,dum2)); 

%  Check  for  convergence 

if  max(max(isnan(w)))==  1  ;error('  it  blowed  up ! ! ! !  ');end 
%  Save  the  variables  to  a  file 
save  tez.mat  c  w  Re  ops  r_inf  cf  kpi  L  N  M  dt; 
end  %  end  of  time  loop 
%%%%% 

function  [Dl,D2,D3]=D_Bmat(M,N) 
global  M3  R  dt  Re 
k=(M+l)*N; 
Dl=zeros(k,M+l); 
D2=zeros(k,M+l); 
D3=zeros(k,M+l); 
fori=l:N 
forj=i:N 
ifi==j 
Dl((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(- 1  *(i- 1  )*(i).  *eye(M+ 1  )-Re/dt.*R.  A2); 

D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(-l*(i-l)*(i).*eye(M+l)+Re/dt.*R.A2); 

D3((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 
=M3+(-l*(i-l)*(i).*eye(M+l)); 
end 
end 
end 

%%%%%%% 
function  E=emat(M) 
E=zeros(M,M); 
m=0:M; 

z=cos(pi.*m./M); 
for  i=0:M 
forj=0:M 
ifi~=j 
if  ((i==0  I  i=M)  &  G==0  I  j==M))  I  ((i~=0  &  i~=M)  &  (j~=0  &  j~=M)) 

c=l;else 
if  (i==0  I  i=M)  &  (j~=0  &  j~=M) 
c=2;else 
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if  (i~=0  &  i~=M)  &  (j=0  I  j==M) 
c=l/2;end;end;end; 

E(i+lj+l)=c*(-l)A(i+j)/(z(i+l)-z(j+l)); 

else; 
ifi==0 

E(i+l,j+l)=(2*MA2+l)/6; 
else 
ifi=M 

E(i+l,j+l)=(2*MA2+l)/(-6); 
else 

E(i+l,j+l)=-l*zG+l)/(2*(l-z(j+l)A2)); 
end 
end 
end 
end 
end 
%%%%%%% 

function  R2=nonlin(w,c) 

global  kmat  rsquare  cotteta  N  M  Ml  c_  c_teta  c_r  teta 

%  variable  list 

% 

%  dw/dteta=nonl  dc/dteta=non2  dw/dr=non3  dc/dr=non4  (in  phy) 

% 

%  transpose  matrices  for  column  operations  in  teta  direction 

r=sqrt(rsquare)'; 

w=w';c=c'; 

nonl=operat(w); 

non2=operat(c); 

wph=phy(w);cph=phy(c); 

non3=[Ml*wph']'; 

non4=[Ml*cph']';non4(:,M+l)=-0.5.*sin(teta); 

term  1  =- 1  ./r .  *  non3 .  *  (non2+c_teta) ; 

term2=  1  ./r.  *  non  1 .  *  (non4+c_r) ; 

term3=- 1  .*cotteta./r.*wph.*(non4+c_r); 

term4=- 1  .*cotteta./r.*non3.*(c_+cph); 

term5= wph./rsquare' .  *  (non2+c_teta) ; 

term6=(c_+cph)./rsquare'.*non  1 ; 

termO=term  1  +term2+term3+term4+term5+term6; 

ter=[termO;zeros(l,M+l);-l.*flipud(termO(2:N,:))l; 

RR2=fft(ter); 

R2=imag(RR2(l  :N,:));R2=R2'; 
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%%%%%%%%% 

function  y=operat(x) 

global  N  M  kmat 

re=zeros(N,M+l); 

X=re+i.*x; 

Y=i.*kmat.*X; 

Y_ful=[Y;zeros(l,M+l);flipud(Y(2:N,:))]; 

y=ifft(Y_ful); 

y=real(y(l:N,:)); 

%%%%%%%%% 

function  y=phy(x) 

global  N  M 

re=zeros(size(x)); 

X=re+i.*x; 

Y_ful=[X;zeros(l,M+l);conj(flipud(X(2:N,:)))]; 

dum=ifft(Y_ful); 

y=real(dum(l  :N,:)); 

%%%%%%% 

function  y=solsvdl(u,s,v,b) 
w=diag(s); 

wmin=max(w)*  le- 1 2; 
fori=l:length(w) 

if  w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 
end 

www=diag(ww); 
y=v*www*(u'*b); 

%%%%%%% 

function  y=solsvd(u,s,v,b) 

w=diag(s); 

%wmin=max(w)*  le-8; 

for  i=l:length(w) 

%  if  w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 

ww(i)=l/w(i); 
end 

www=diag(ww); 
y=v*www*(u'*b); 
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